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^ ! Abstract 

g 

This paper describes a novel numerical approach to find the statistics of the 

d ■ non-stationary response of scalar non-linear systems excited by Levy white noises. 

c/) ! The proposed numerical procedure relies on the introduction of an integral trans- 

form of Wiener-Hopf type into the equation governing the characteristic function. 
Once this equation is rewritten as partial integro-differential equation, it is then 
solved by applying the method of convolution quadrature originally proposed by 
Lubich, here extended to deal with this particular integral transform. The pro- 
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Q ■ posed approach is relevant for two reasons: 1) Statistics of systems with several 

O . different drift terms can be handled in an efficient way, independently from the 

kind of white noise; 2) The particular form of Wiener-Hopf integral transform 
CN ■ and its numerical evaluation, both introduced in this study, are generalizations of 

fractional integro-differential operators of potential type and Griinwald-Letnikov 



> 

o 



O , fractional derivatives, respectively. 

cn 

i> ■ 1 Introduction 

o 

Non-linear stochastic differential equations (SDEs) model many real phenomena. Since the 
early works of Brown and Langevin, SDEs driven by Gaussian white noise processes have been 
. studied and for long time they have represented the standard approach in stochastic modeling. 

rN I Only in the last decades. Levy white noise excited systems have attracted the interest of 

C^ ■ many authors in finance, in physics and in engineering. The reason of such an interest is that 

many physical phenomena exhibit non-Gaussian behavior which cannot be neglected: skewed 
probability distributions, inverse power-law decay of the density and divergent moments are 
typical characteristics which cannot be modeled by Gaussian distribution. 

The generalized form of the central limit theorem given by Paul Levy provides the nec- 
essary mathematical support to justify the choice of non-Gaussian excitation, stating the 
conditions under which the sum of independent random variables converges to the so-called 
Levy a-stable distribution. This is a class of distributions defined in terms of the stability 
index < a < 2, being a = 2 the Gaussian distribution. 



*" E-mail: giuliocottone@unipa.it; giulio.cottone@tum.de 



If on the one hand non-Gaussian assumption is useful to enriching models which become 
more and more realistic, on the other hand this increases the mathematical complexity of 
the governing equations. For example, the Fokker-Planck equation, regarding the time evo- 
lution of the probability density function, is a partial differential, an integro-differential or 
a fractional differential equation in the cases of Gaussian, Poisson and a-stable white noise 
external excitation, respectively. Very few analytical solutions are available and numerical 
methods to find the statistics of the solution of the non-linear SDE are frequently the only 
valid alternative. 

In this context, this paper introduces a relevant numerical method for handling with 
this problem in case of externally driven scalar SDEs. It will be shown that the proposed 
method applies to every kind of Levy white noise external excitation and that it is numerically 
advantageous. To frame this method, it is useful to briefly recall some approach given in 
literature, considering separately the Gaussian, the Poisson and the a-stable Levy case. 

SDEs excited by Gaussian white noise have been studied for a long time and, for some 
particular form of the drift and diffusion coefficients, stationary density [201 EZ] are avail- 
able by solving the associated Fokker-Planck-Kolmogorov (FPK) equation. Non-stationary 
solution can also be, at least theoretically, found [32j if the eigenfunctions of the non-linear 
FPK are available or can be computed in numerical way. Among the proposed strategies to 
approach the case of Gaussian excitation we mention: 1) methods to solve the FPK equation 
[211 EHl EH SOI m] by numerical approaches, in particular by Finite Element approach; 2) 
path integral method [2H]; 3) methods to transform the original non-linear systems into an 
equivalent systems, which allow the approximation of high dimensional systems, such as the 
stochastic hnearization f33] and the stochastic averaging fTS]. While methods 1) and 2) are 
confined to dynamical systems with relatively small number of state variables and are still 
impracticable [2S] in high dimensions, methods in 3) may not be accurate in case of strong 
non-linear drift. 

For Poisson white noise excitation, generalizations of the Fokker-Planck equation for the 
evolution of the response density, indicated as Kolmogorov- Feller (KF) equation, are given in 
[HI [17] and they are based on a generalized form of the Ito's calculus; another approach 
starting from the theory on filtered Poisson processes [13] leads to the equations governing 
the characteristic function in case of polynomial drift and external Poisson white noise, that 
are the spectral counterpart of the KF equation. Exact stationary solutions for few cases have 
been proposed in [30j. Numerical procedures to find statistics [Ij, equivalent linearization 
[12], path integral solution [9l [29] and extensions of the Bogolyubov's theorem (stochastic 
averaging) [19j are also available. 

The study of Levy white noise processes is characterized by mathematical complexity 
mainly due to the heavy tails of the a-stable distribution, which imply that moments of order 
p > a are divergent. This causes the inapplicability of many standard procedures of stochastic 
calculus, based on moments. Ito's calculus has been extended to deal with Levy white noise 
processes by many authors and readers are referred to ^31j for a rigorous mathematical treat- 
ment and to [15] for a complete presentation of the topic with applications in sight. Although 
this topic attracted interest only in the last decades, a considerable amount of research has 
been conducted: analytical stationary solution for particular polynomial potential are given in 
[21 [3 [ini [TT] ; equivalent linearization technique based on the minimization of a functional of 
the characteristic function is presented in [l3] ; methods based on the continuous time random 
walk leading to the generalization of the Fokker-Planck equation in terms of fractional Riesz 
derivatives can be found in the topic reviews [251 EI]- 

To the author's knowledge none of the above cited method is applicable, without substan- 
tial modification, to the treatment of SDEs excited by any external Levy noise. This paper 
fills this gap showing a numerical method to find the statistics of the non-stationary response 
of non-linear systems driven either by Gaussian, Poisson or a-stable white noise processes. 



The solution strategy is based on a proper numerical treatment of the equations ruling 
the characteristic function that are indicated as Spectral Fokker-Planck equations (SFPEs). 
Main advantage of these SFPEs is that their mathematical form remains unchanged depend- 
ing on the excitation. Although dealing with the characteristic function rather than with the 
probability density is a strategy already used by other authors, see [3J, an unique numerical 
treatment for any Levy white noise has never been proposed, to the best of the author's knowl- 
edge. Moreover, it is important to stress that this method is not restricted to drift coefficient 
of polynomial type, or to power-law type as in [3], but conversely is extremely advantageous 
and straightforwardly applicable to any non-linear form of the drift, as the numerical section 
shows. To this aim, we will give firstly a computationally suitable form to the governing equa- 
tions in order to highlight the integro-differential structure of the equations. Then, once the 
integro-differential nature of the SFPEs is made explicit, a Wiener-Hopf integral transform 
will be introduced and the convolution quadrature method, originally proposed by Lubich 
[211 [221 ESI, is reformulated in order to handle such an operator. In this way, the character- 
istic function of the response of non-linear stochastic differential equation is found solving a 
simple linear system of ordinary differential equations. 

It is worth to mention that one of the strength of the proposed approach relies on the 
use of convolution quadrature method, that here is properly modified to suit our problem in 
a straightforward way. Indeed, as pointed out in [36] convolution quadrature has excellent 
stability properties compared to other discretization method and is a valid alternative both 
from the accuracy and the efficiency point of view. 

2 Problem statement 

Consider the nonlinear system excited by an external white noise process, in the form 

r x{t)^f{x{t),t) + WL{t) 

^' \X{0) = Xo 

where / (X (t) , t) is non-linear deterministic function of the response process X{t) and Xq is an 
initial condition, which might be either deterministic or random, with assigned distribution. In 
the framework of the generalized theory of random processes a stationary white noise process 
Wl(^) is defined as the formal derivative of a process with stationary orthogonal increments 
L{t) and starting from zero a.s., i.e. L(0) = 0, called Levy process. Formally we can write 

that, introduced into ([T]) gives the Ito's form 

(2) dX (t) = / (X (t) ,t)dt + dL (t) 

In what follows we restrict our attention to the following three interesting cases of external 
excitation: 

1. Gaussian or normal white noise process Wait) defined as the formal derivative of the 
Brownian motion B(t) whose increments dB{t) — B{t + dt) — B{t) have zero mean 
Gaussian distribution and with constant power spectral density equal to q/2Ti^ where q 
indicates the noise strength, assumed constant in time. 

2. Poisson white noise process Wp{t) defined as the formal derivative of the Compound 
motion C(t); Wp{t) is defined as sum of Dirac's delta impulses S{t — Tj^), having random 



distributed amplitudes Y/^, occurring at Poisson distributed random times Tj^ indepen- 
dent from Yfc, explicitly 

N(t) 

k=i 

where N{t) indicates a counting process with mean arrival rate A. In the following, 
it is assumed that the distribution of the random variable Y is known along with its 
characteristic function (J)y{0). No further assumptions are done regarding the existence 
of moments of Y . 

3. a-stable Levy white noise process Wa{t) defined as the formal derivative of the a-stable 
Levy motion La{t)^ see [15]. 

By characterizing the driving noise, it is possible to derive the equation governing the 
evolution of the characteristic function, that we indicate as Spectral- Fokker-Planck equations 
(SFPEs) (see [5j). For stochastic differential equations excited by external Gaussian, Poisson 
and a-stable Levy white noise processes, SFPEs are written, omitting arguments, as 



(3a) 



iOE 



f{X,t)e 



iex 



(3b) 
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f{X,t)e 



iex 



A0(l-0y) 



(3c) 



i9E 



f{X,t)e 



idX 



\o\' 



where i = ^/^ is the imaginary unit, E [•] denotes expectation and (/) = E [e^^^^^^j is the 
characteristic function of the response X{t). Readers are referred to [TSl ES] and to [8J for 
different prospectives to derive the characteristic function equations of the state of dynamical 
systems driven by Levy white noises. 

From the above SFPEs it is noticed that the non-linear drift term influences the time 
evolution of the characteristic function only in the term E [f (X, t) e^^^^^^j . 

The case in which / (X, t) is polynomial of the state variable X(t) is handled easily. In 
such a case Eqs.(j3]) become partial differential equations by the property 



(4) 



09^ 



(iX) 



k ^iOX{t) 



which follows from the definition of (j){9^t) and standard numerical schemes for the solution 
of partial differential equations can be used. 

On the contrary, if /(X, t) is not a polynomial in X(t), the term E [/(X, t)e^^^*^^)] is not 
directly related to the function (j){0^t). 

Aim of this paper is to provide a general numerical framework to solve the Ito's stochastic 
differential equation ([2]) in which the non-linear drift term is not of polynomial form. The 
strategy used is based on the solution of eqs.(|3]). First, the expectation E [/(X, t)e^^^(^)] is 
rewritten as a Wiener-Hopf convolution integral operator applied to the characteristic function. 
Then, the Lubich's method for the quadrature of convolution integral will be adapted to find 
a form that is numerically convenient. Finally, the general problem of numerical evaluation of 
such an integral by convolution quadrature will be addressed and a procedure for the numerical 
evaluation of Eqs.® will be proposed and tested by some examples. 



3 Convolution quadrature for the integro-differential 
equation governing the characteristic function 

Indicate with X{t) the response to the stochastic differential equation in ([T]) and let p{x) 
and (j){9) denote the density and the characteristic function of X(t), respectively, where the 
dependence on time is for readability's sake not explicitly indicated. 

From the definition of the expectation operator and recalling that the characteristic func- 
tion is the Fourier transform of the density p(x), the expectation E [f (X) e^^^] can be written 
as 

J —OO 
-| POO /♦OO 

= — / / </>(«) e-'^'^du f (x) e'^'^dx = 

^'^ J — OO J — OO 

-1 /'OO roo 

(5) =— / cf>{u) f{x)e-'^''-^^''dxdu 

^'^ J —OO J — OO 

The inner integral is the Fourier inverse transform of /(x), which can be indicated by f{u)^ 
i.e. 

(6) /(«) = i_|^/(^)e— dx 
Introducing the latter in eq.([5]), the relation 

(7) E f{X)e'''' = / <Piu)fiu~9)du 

«^— OO 

in terms of the characteristic function 0(0) is obtained. This is an integral transform that 
can be referred to as Wiener-Hopf type, in analogy to the Wiener-Hopf integral equations. 
For what follows, it is more convenient to dub the previous convolution integral as (Tfcj)) (6) 
that is an integral transform applied to the characteristic function (f){9) and depending on the 
function f{x). 

The advantage of this representation relies on the easy and efficient numerical evaluation 
of the transformation operator {Tf(j)){9) by the convolution quadrature shown in the following. 
This is adapted from the original work of Lubich [2T1 [22l |23], which readers are referred to. 

Let us indicate by 



and 



roo 

F-{s)^ / f{-t)e-''dt 
Jo 

roo 

F+(.)= / f{t)e-^'dt 
Jo 



the Laplace transform of the function f(t) for negative and positive values of the variable t, 
respectively. Moreover, let us recall the inverse Laplace transform relations 

(8a) f^t) = ^J^F+{s)e''ds 

(8b) f^-t)^^.j F.{s)e''ds 



The latter can be introduced into the definition of the operator {Tf(f)){9)^ and after some 
algebraic manipulation one obtains 

(9) [Tfci)){e)^ f^-t)<i){e-t)dt+ f{t)<p{e + t)dt 

Jo Jo 

Introducing (|8j) in the latter and rearranging, produces 



(10) {Tf^) {9) = ^f F_ (s) f {u) e^^-^^^duds^ 

— oo 

oo 

+ — I F^{s) f c^{u)e<''-^^duds 
2tii Jy J 

e 

The inner integral in the first term can be thought as the formal solution to the differential 
equation 

y'{e)^-sy{e)-<p{e) 



^^^) \y(-oo) = 

whose solution will be denoted by y{9^ s) to remark the dependency on the variable s of the 
Laplace transform. In similar fashion, the second term can be expressed as the formal solution 
to the differential equation 

z'{e)^-sz{e)-(j){e) 



^^^^ I z(oc) = 

whose solution is z{9^s). Eq. (lTT]) can be rewritten as 

(13) {Tf<p) (^) = ^ / F- {s) y (9, s)ds + ^.j F+ {s) z {6, s) ds 

It is worth to note that the term E [/ {X) e^^^] , taking into account of the non-linearity 
of the system, is thus expressed as an inverse Laplace transform of a functions depending of 
the unknown characteristic function (j){9). The main advantage of this representation relies on 
the suitable numerical solution that can be performed following the work of Lubich [21} [22] 
on quadrature of convolution integrals. The main point, outlined in the following, is that the 
solutions to the differential equations in eq. (ITTI) and eq. (fT2]) are expressed in terms of sums by 
applying difference schemes. The following way to derive a series approximation slightly differs 
from the original method of Lubich that uses formal power series and operational calculus. 
Main motivation is that eqs. dlip and (I12p have an initial condition at infinity that can be 
taken into account by using recursive relations as shown below. 

Eqs. dlip and (I12p will be approximated by means of a simple Eulerian scheme for the first 
order derivative. Higher approximation might be pursued by higher order derivative schemes, 
but to keep a simpler notation, we will confine our attention to the simple Eulerian one. 

Let us indicate 9^ — nh, y {9^) — Vm ^ i^n) — ^m (^n) — <t^ni where /i G M is a grid step 
and n G N. 

The approximation of eq. (fTT]) can be pursued by a forward Eulerian scheme that reads 

(14) Vn - Vn-l = Shyn + (pnh 

To find a solution to the latter equation, define the infinite dimensional vector Y = 
[ yn yn-i ••• ] 5 collecting the solution at sampling steps, the infinite dimensional vector 
^ — [ 4^n 0n-i ••• ] Q^nd the formal matrix 

1 
1 




such that the infinite set of equations represented by eg. (1141) can be written in the form 

(15) Y = (1 - shy^ AY + */i (1 - sh)-^ 
Formal solution to the latter is provided in the form 

(16) Y = Tl - (1 - sh)-^ a) "^ */i (1 - sh)-^ 

being 

" 1 ••• 
1 



an infinite dimensional identity matrix. Because of the particular structure of the matrix 
(l — (1 — sh)~ A), it is not hard to prove by induction that its inverse has the form 



(17) 



(l - {1 - sh)-^ a) 



1 {1-shy^ {l-sh)~^ {l-sh)~^ 
1 {l-sh)~^ {l-sh)~^ 

1 (1-s/i)"^ 



Concluding the reasoning, the general term of the differential equation yn follows from 
eqs. (fT6]) and (IT7|) and after some algebraic manipulation it can be written as 



(18) 



vn = E(-i; 



k+1 



k=0 



(t^n-kh ^ 



In the same way, we can handle eq. (fT2l) , but, in this case, by applying more conve- 
niently a backward Euler scheme because of the different boundary condition. Calling Z = 
[ Zn Zn+i ... ] , $ = [ 0ri 0n+i ••• ] it is casy to show that the resolvent equation can 
be formally expressed in matrix notation as 

(19) Z = Tl - (1 - sh)-^ a) "^ */i (1 - sh)-^ 



and that the general term Zn is found to be 

oo 

(20) zn = j2 (-i; 



k+1 



k=0 






Introducing the expression found for y^ and Zn into eq. (fT3l) and performing the integral 
inside the sum, one has to evaluate the two integrals 



(21a) 



-I 

27TI Jy 



F-{s) 



r (5 - 1/h) 



k+l 



1 d^F_(s) 
as = -— 



k\ ds^ 



OikU) 



s=l/h 



(21b) 



27rz Jy 



F+(s) 



r{s- 1/h) 



k+l 



& = i'''^+W 



k\ ds^ 



s=l/h 



^k{f) 



due to the Cauchy's integral formula for the derivatives. The notation adopted to indicate 
the coefficients ak(f) and Ok{f) is meant to underline their dependence on the function /. 
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Sumining up, from the non-linear function f{x) we obtain / performing the inverse Fourier 
transform and by eqs.(j8j) the functions F±(s) are calculated. These functions are expressed 
by Taylor expansion around s = 1/h and the coefficients are ak ^tnd cjk- Once the coefficients 
are known, (11 3p is finally expressed in the form 

(OO CXD \ 

J2Mf)^{0 - kh) + 5^a;,(/)0(0 + kh) 

having indicated with ak{f) = (-1)^+^ h'^akif) and with ojk{f) = (-1)^+^ h'^Ukif). 

If one selects a fixed small value of the increment h in the latter expression, it corresponds 
to a first order approximation of the operator {Tf(j)){9)^ otherwise, the expression is exact to 
the limit. 

3.1 Examples 

Some simple examples will show how to practically apply the relations found. In this section, 
for readability's sake, the dependence of the coefficients on the function will not be indicated 
because evident from the contest. 

3.1.1 SDE with linear drift. 

The representation formula of the non-linear drift given in eq. (l22]) is valid under the hypothesis 
that the involved integral transforms exist, at least in the generalized sense. Let us first 
consider a simple linear system excited by external Gaussian noise, that is with a{X) — a — 
constant and h{X) = 1 

a > 



(23) dX (t) = 


-aX (t) + dB (t) 




The SFPE equation in this case is 






t 


= -iaOE 


'Xe''''' 


^2 

2 



which becomes a partial differential equation once is recognized that E [Xe^^^] = ~^af ' ^^ 
stated in eq.(jl|). We follow the procedure proposed to attain to this simple result. 

In this case, f{x) = —ax, and the inverse Fourier transform is f{0) = —ia5'{9), where 5 is 
the Dirac's delta and the prime indicates the first derivative, and F±{s) — ^ias/2 by direct 
integration. The integrals in the Laplace transforms have been performed with lower limit in 
0, and not in 0—, and for this reason the factor 1/2 appears. It is easy to find the Taylor 
series in 5 = 1/h for both functions F±{s) 

ia ia ( 1 

+ ^ ^ 2h 2 V h 

_ , , ia ia f 1 

and consequently to calculate the coefficients that are 

ia _ ia _ _ 

ao^-—-; «! = ---; ^2 = Q^3 = ••• = 

_ ia _ ia _ _ 

^0 = ^5 ^1 = y ; CJ2 = CJ3 = ... = 



Lastly, introducing the latter coefficient in eq. ([22|) and making the limit h ^ we obtain 



■'"w 



that is the searched result. Without taking the hmit, the construction of the operator {Tf(j)){9) 
produces the first order approximation of d(j)/d9 by finite differences. 

3.1.2 SDE with non-linear polynomial drift. 

Generahzation to non-hnear drift term of power-law type, i.e. f{x) — x^ , with j G N follows 
plainly. Indeed, in this case f{9) = i^ 5^^\9)^ where the apex between brackets indicates the 
jth derivative and F±{s) — {±iy s^ . Substitution in eq. (l21ap and (I21bp gives 

that, introduced in ([22]) gives {Tf(l)){9) = (-i)^ c/)^^^ (6>), for h ^ 0. 

3.1.3 SDE with real power-law non-linear drift. 

Consider a drift term of the form / (x) = \x\^ sgn (x). Non-linearity of this kind is frequently 
used to model viscous-elasticity. Levy driven non-linear system with this kind of drift have 
been studied for example in ^ and ^ by applying the fractional calculus. Application of the 
procedure proposed leads to interesting results, in accordance with the above cited papers. 
First of all, by means of Mathematica, it is found that 

zcos (77r/2) r (1 + 7) |^i_i_^ _ /z)\ 7:^ / \ I ^^^ 



The coefficients 

«. - .r-^^h-r-^ (71; -. - - . ■ ; ,.. h^-' ' ' 



2sin(77r/2) V ^ / ' 2sin(77r/2) V ^ 

introduced in f[22l) give that the integral {Tf(j)){9) coincides in this case with a fractional 
operator of potential type (see [34j, p. 214), that is 

iTfcf>m ^ (if->) (6) = ^— ^ {{Dlcf>) (9) - {DU) {6)) 

where {H~^(j)) (9) is called complementary Riesz fractional derivative of real order 7 > 0, 
7 7^ 2, 4, .... It can also be defined in terms of Griinwald-Letnikov fractional derivatives in the 
form 

(25) {Dld^) {9) = hmh-rf^i-l)' ( ^ ) ^PiO^kh) 

and such coefficient coincides with those found by the proposed convolution quadrature. 

Remark: The previous examples show the relations between the transform operator 
{Tf(j)){9) and the classical derivatives and the fractional derivatives. 

The first two examples showed that {Tf(j)){9) coincides with the classical derivatives of 
integer order fc, with A: G Z in case the function /(x) = x^ is selected. 

If we consider / {x) — \x\~^ or / (x) = |x|~^sgn(x), with 7 G C, then the third example 
shows that the integral transform {Tf(j)){9) coincides with a particular Riesz fractional integral. 
Moreover, the numerical procedure based on the convolution quadrature through the Eulerian 
schemes to solve (fTTI) and (fT2l) gives: 1) the finite difference schemes of first order, in case 
of integer power-law kernel; 2) the Griinwald-Letnikov fractional approximation in case of 
complex power-law kernel. 

This similarity suggests that many other properties typical of fractional calculus, such 
as the integration by part formula, the composition rules and the definition on a bounded 
interval, may be stated also for the integral (Tf(j))(9). We report such properties in Appendix 
A, for readability's sake. Conditions for the convergence of the sum are reported in Appendix 
B. 

9 
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4 Solution of scalar non-linear stochastic systems 

Let us revert to the statistic of the solution of the non-hnear systems ([T]). Introducing the 
transform {Tf(j)){9) in (|3]) the integro-differential equations in terms of the characteristic func- 
tion (j){9^t) can be written in the form 

(26) ^^^^ie(Tf<i>){e)-g{6)<l^{e,t) 

in which g{9) is equal to 0^/2, A(l — 0^(0)), l^l*^ in the case of Gaussian, Poisson and a-stable 
external white noise, respectively. To solve this equation, let us select a finite interval [— ^, 6] 
and divide it in equidistant subinterval of finite dimension h — 9/m^ with tti G N. The choice 
of a symmetric interval depends on the symmetry of the characteristic function. Moreover, 
indicate by 9j — jh, (j)j — (j){9j)^ with j — —m,...,m. Formula to express {Tf(j)){9) in the 
bounded domain [—9^9] is given in Appendix, eq. ([37l) . Then ([26|) is rewritten as 

/[^] [^] ^ 

k=0 k=0 

V / 

As the latter equation is valid for every j = —m, ...,7n, the linear system of differential 
equations with the form 

(28) = R(/) 

can be constructed, in which R = iUTf — G, 

(29) U= r-0^,...,O,...,0^J 
is a diagonal matrix, 

/ ao (/) + C^O (/) OJi (/) U2 (/) U2m (/) \ 

C^l (/) C^O (/) + ^0 (/) ^1 (/) • • • ^2m-l (/) 

(30)Tf= a2{f) aiif) cq (/) + cc;o (/) '•. 

\ Oi2m (/) Oi2m-l (/) <^1 (/) <^0 (/) + ^0 (/) / 



^ - [0-m,...,0m] 



is a Toeplitz matrix of coefficients depending on the non-linearity, and 4> 

is the vector collecting the 2m + 1 unknown values of the characteristic function at the grid 

points. 

The matrix G is 

1. G = 1/2 1"—^^, ..., 0, ..., 0^J for Gaussian noise; 

2. G = A [1 — (/>y (— ^m), ..., 0, ..., 1 — 0y (0m)J foi" Poisson noise; 

3. G = A [|— ^mT 5 •••5 O5 •••5 l^mTJ foi" Q^-stable Levy noise. 

In the stationary case (/) = 0, and the system of differential equations (|28]) becomes 
an algebraic linear system, that can be solved, by posing the relevant boundary condition 
00 = 0(0) = 1. Such a boundary condition can be easily posed by defining the reduced 
2m dimensional vector of coefficients b-^ = [ ujm ••• ^1 <^i ••• <^m ] and the reduced 
matrix R obtained by deleting both the central line and the central row to the matrix R. In 
this way the reduced vector of unknown 4>^ — [0-m, •••5 0-1, 0i---0m] is the solution of the 
simple linear system in the form 
(31) 4) = -R-^ b 

The non-stationary solution can be found by solving the linear system of differential equa- 
tions (|28]) . This can be simply achieved by spectral decomposition. 
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5 Validation and Applications 

To assess the accuracy of the proposed procedure we start studying the SDE in the Ito's form 

(32) dX (t) = {-X (t) + cos {X (t))) dt + dB (t) 

in which dB{t) is the standard Brownian motion, with q — 1. In the stationary case, the 
statistic of the solution are known in analytical form and they can be used to validate the 
proposed method. In fact the stationary density is given as 



Pi 



(^x^^Ce^foi-i+^osm 



being C a normalization constant (see [37j). 

For this application we selected ^ = 15, and different mesh sizes m = 100; 150; 200; 500 
to investigate on the convergence of the numerical approximation to the exact solution and 
the influence of h. The linear algebraic system of equations given in (|3T|) , which returns 
the stationary characteristic function, has been solved and then by inverse numerical inverse 
Fourier transform, the stationary density has been found for every value of h and it is indicated 

^8 Ph{x). 

The coefficients defining the numerical form of the integral transform {Tf(j))(9) are 



and 



2 I r(i + j) 







2 \ r(i + j) h 



where Tj = 1 for j = and Tj — h ioi j ^ 0. Figure ([T]) shows the absolute error, calculated 
in the positive axis, between the exact and the approximate stationary density of the solution 
Ph(x). This result shows that by this method it is possible to achieve a good numerical 
approximation keeping an absolute error that is at most of order h. This order of magnitude 
derives from having used an Eulerian scheme in finding the series approximation of the integral 
transform {Tf(j)){9). In order to reduce such an error, higher order approximations can be used 
and their implementation is underway by the author. 

The non-stationary density of this system has also been calculated for m — 100 by solving 
eq. (l28]) . For this case of external Gaussian excitation, the solution of the Fokker-Planck 
equation can be approximated also by eigenfunction expansion (see for instance [32J, p. 101). 
In Figure ([2]) a comparison in term of first and second order moments calculated from the 
non-stationary densities calculated by the proposed method and the eigenfunction expansion 
method are compared, showing good agreement. 

After having validated the method, in the following, we consider the non-stationary solu- 
tion of the three non-linear SDEs: 

dX it) = {-X (t) + cos {X (t))) dt + dL (t) 
X (0) = Xo 

in which: 

1) for dL{t) = dB{t), the initial condition is Gaussian, i.e. Xq = N (3, 1) and m = 100; 

2) for dL{t) = dC{t) we consider exponential distributed jumps with mean /xy = 3, Ay = 3, 
Gaussian distributed initial condition Xq = N (0, 1) and m = 100; 

3) for dL{t) = dLa{t), i.e. in case of a-stable Levy excitation, we selected a = 3/2, uniform 
initial condition Xq = t/(0, 1) and m = 100. 
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10.0 



Figure 1: Absolute error for different values of h: h = 0.15 circles, h = 0.10 squares, 
h = 0.075 rhomboidal and h = 0.030 triangles 




(a) 



(b) 



Figure 2: Time evolution of the first and second order moments. Comparison between 
the results calculated from the eigenfunction expansion of the FPK (continuous line) 
and the proposed method (dotted) 



It is usefull to recall that the Fokker-Planck equations, governing the densities, for these 
systems have the forms 



(33a) 
(33b) 

(33c) 



dp 
~dt 



d_ 
dx 



dp d 1 dp 

/OO 
p{Opy {x - ^) 
-OO 



-Od^ 



dp 
~dt 



dx 



({-X + cos (x)) p) + (P» (x) 



in which V^ is the Riesz fractional derivative, see [26j for the derivation of the equation and 
[M] for relevant properties of this operator. The numerical solution of the latter equations is 
not trivial: In the case of Gaussian excitation, in which the eq. fl33ap is a partial differential 
equation, numerical evaluation of the eigenfunctions and eigenvalues of the Fokker-Planck 
operator for obtaining the non-stationary solution is required (see [32], ch. 5); Eqs. (l33bp and 
fl33cp are respectively a partial integro-differential and a partial fractional differential equations 
and quadrature methods can be applied. 

The SFPK equations, that are their spectral counterpart, are instead partial integro- 
differential equations of the same kind 



(34a) 



i9E 
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(c) 



(d) 



Figure 3: Gaussian white noise excitation: a) Real part of the CF; b) Imaginary part 
of the CF; c) Contour plot of Re{(f){9,t))] d) Contour plot of Im{(f){6,t)) 



(34b) 



^iOE 



ex 



(-X + cos(X))e^^^ -A(/)(l-(/)y) 



(34c) 



i9E 



(-X + cos(X))e^^^ ] 



and, by the procedure previously shown, are simply rewritten as 

(35) (^ = iUTf (/) - Gcj) 

The system of linear differential equations in eq. ([35l) under the approximated boundary con- 
ditions (j){9^ t) — and proper initial conditions, thus obtaining the approximation of the time 
evolution of the solutions' characteristic functions for each of the non-linear system. Of course, 
the coefficients defining the matrix UTf are calculated only once because dependent only on 
the drift. These characteristic functionals are complex functions and Figures ([3]), (jlj) and ([5]), 
show both the real and the imaginary part. 

It is worth noting that eqs. (f3T]) might be solved with other finite difference rules or quadra- 
ture formula, but the Lubich's convolution quadrature, core of the method here provided 
exhibits in general higher numerical stability, accuracy and efficiency 
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(a) 



(b) 




10 



-10 





(c) 



(d) 



Figure 4: Poisson white noise excitation: a) Real part of the CF; b) Imaginary part of 
the CF; c) Contour plot of Re{(j){e,t))] d) Contour plot of Im{(j){e,t)) 

6 Conclusions 

This paper has shown a numerical procedure to find the statistics of non-linear dynamical 
systems under Levy white noise processes. 

First, the Spectral Fokker-Plank equation, which governs the time evolution of the charac- 
teristic function, has been rewritten in terms of a Wiener-Hopf integral, indicated as {Tf(j)){9). 
This term depends only on the non-linear drift. For every white external excitation, either 
Gaussian or Poisson or Levy a-stable, the mathematical structure of the governing equations 
in terms of the characteristic function remains the same. 

This method exploits this fact, providing an efficient scheme to evaluate the Wiener-Hopf 
integral: to this aim the Lubich's convolution quadrature has been reformulated to deal with 
this particular case. The partial integro-differential equations of the characteristic function 
have thus been transformed into a linear system of ordinary differential equation that can be 
easily solved. In this way the non-stationary solutions of three non-linear stochastic systems 
have been found and presented to support the method. 

Validation and error analysis has also been presented by comparison with a benchmark 
solution in case of external Gaussian excitation. Moreover, the convergence of the numerical 
scheme on which the method relies on has also been provided in Appendix B. 

Summing up, the advantages of the proposed method relies on: i) excellent numerical 
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(a) 



(b) 



9 0.0 




0.0 




(c) 



(d) 



Figure 5: a-stable Levy white noise excitation: a) Real part of the CF; b) Imaginary 
part of the CF; c) Contour plot of Re((j)(9^t))] d) Contour plot of Im{(j){6^t)) 

stability given by the Lubich's convolution quadrature method; ii) applicability to every kind 
of external white noise excitation; iii) ease in constructing the matrices defining the system of 
linear differential equation that have Toeplitz form. 

From the theoretical point of view, it is worth to make some further considerations on the 
form of the Wiener-Hopf integral {Tf(j)){9) and on its series form found in (j22l) , in connection 
to fractional differential operators. 

For a long time fractional derivatives and integrals have been considered just particular 
integral transforms with the property of interpolating ordinary derivatives. Properties with 
respect to Fourier transform, the Griinwald-Letnikov form, the compositions rules, and many 
others have then supported the use of such a mathematical tool that nowadays is widespread. 

This paper deals with the integral transform {Tf(j)){9) in perfect analogy with the frac- 
tional integrals, introducing similar notation and giving the most important properties. This 
simplifies the use and the mathematical treatment of the specific problem at hand on SDEs. 
In appendix A the similarities between the most useful properties of fractional calculus and 
the correspondent properties for the Wiener-Hopf integral transform have also presented. 
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A Considerations on the integral transform (Tf(j)){6) 

The integral transform {Tf(j)){9) defined in ([7]) is a particular type of convolution integral in 
which the kernel is the inverse Fourier transform of a given function /. Due to this particular 
structure of the kernel, many similarities between this integral transform and fractional inte- 
grals can be drawn and can be useful in applications. In this appendix, by simple application 
of standard theorems on Fourier transformable functions, these properties are reported. 
First, we have shown it this paper that the integral transform 

/oo 
<t>{u)f{u-e)du 
-oo 

by means of a modification of the Lubich's approach can be expressed as 

(oo oo \ 

y^ak{f)(t>{0-kh) + y^uk{f)(t>{e + kh)\ 
k=0 k=0 J 

which reminds the relation between the Riemann-Liouville fractional derivatives and their 
series approximation by the Griinwald-Letnikov's formula. 

{Tf(j)){9) has been defined as integral in the whole real domain, but as like as the right 
and the left side Riemann-Liouville fractional integrals, it can be also defined in the semi- 
axis. Indeed it might be useful to split {Tf(j)){9) in two operators in which, the first takes 
into account of all the values preceding 9, indicated by {Tf^(j)){9)^ while the second takes the 
values following 9 and is indicated as {T f _(j)){9). More explicitly, the series approximation of 
these two integrals is plainly obtained by the results of the previous sections and reads 

oo 

(36a) (Tf^c^) {9) = limj^ «,(/)(/) (0 - kh) 

oo 

(36b) (Tf_c^) (9) = ]imJ2^k(m (0 + kh) 

It is plain that (T/^) (9) = (T/^^) (9) + (TfA (9). 

Furthermore, the similitude with the fractional operators suggests how to calculate trun- 
cated versions of the operators (Tf^cf)) (9) and (Tfcf)) {9). Let us consider a bounded in- 
terval [9i^9f]. It is easy to demonstrate that, in perfect analogy with the Griinwald-Letnikov 
fractional integral defined on this bounded interval, the following relations 

[^] 

(37a) (r/^(/.) (9) = lim J^ a^ (/)</. (^ - fc/i) 
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(37b) (r;_0) (0) = lim Yl uk{f)ct>{0 + kh) 

hold true, where the symbol [•] indicates the integer part. 

Other relevant properties can be derived by plain application of the Fubini's theorem. The 
notation introduced for the integral {Tf(j)){9) allows to simplify the following relations. 

If we indicate the Fourier transform as J^ and apply it to (Tf(j))(9) we obtain 

(38a) T{{Tf<t>)ie);u;}^T{cf>ie);u;}f{-u) 

(38b) ^-1 {(TfcP) (9) ; u} = J-^ {</. [Q) ; u} f (u) 

The Fubini's formula allows also to write the equivalent form of the so-called fractional 
integration by parts formula ([34j, p. 96), that in this case assumes the form 

/oo roD 

(TfcP) (e) g (6) dO = / {Tf^g) {9) <P (9) dO 
-OO J —oo 

where /* = /(— x). 

Lastly, analogous form of the composition rule for fractional operators is reported. It can 
be shown that, by applying the integration by parts formula, above reported, and the Fubini's 
theorem, the following formula 

/oo 
cf>{u)Jg{u-9)du 
-OO 

holds true, where 

1 foo 

f9{0)-^J_ g{t)f{t)e-'''dt 

is the inverse Fourier transform of the product of the functions f{t) and g{t). 



B On the convergence of (1221) 



Let us assume that f{x) is Lipschitz and satisfies some growth condition such that the SDKs 
in ([1]) have an unique solution. Moreover, it is assumed that f{x) is Fourier transformable, at 
least in generalized sense, and that the inverse Laplace transforms (|S]) are defined for every 

S > So. 

The only approximation introduced in section 3 is the Eulerian scheme in equation (II 4p 
that in its general form can be written as 

dOVn + Ciiyn-l = hbo(j) (yn) 

and it is known to converge (see, for example [42j, p. 494) if cj) is Lipschitz, i.e 

\(l) (6) - (I) {u)\ <L\e-u\ 

and 

hL\bo\ ^^ 

l^ol 
It suffices to recall that in our case ag = 1 — sh, ai = — 1 and 6o = 1 and to recall that the 
function (f){y) is a characteristic function, i.e. \(j){y)\ < 1 for every real y to derive that the 
sum in (fTSl) converges for 

1-h 

So < s < ——- 
h 

From this relation a proper value of the increment h can be selected for numerical imple- 
mentation. 
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